Seed rain, seed bank, and vegetational dynamics of north-central Missouri remnant and reconstructed tallgrass prairies
Background:
Author: Katherine Wynne (Wynnekat@msu.edu)
Co-authors: E. Eyerly, L. Sullivan
Created: 14 June 2023
Files:
Cleaned_Spring_Cover_Data_July_2019.xlsx - contains spring p/a cover data
Cleaned_Fall_Cover_Data_September_2019.xlsx - contains fall % cover data
Seed_Rain_Data_2019_FINAL.xlsx - contains seed rain data
Seed_Bank_2020_Cleaned.xlsx - contains seed bank data
SR_Traits_Dataset.xlsx - contains trait data (work in progress)
R Version: R 4.3.2 (2023-10-31) – “Eye Holes”
RStudio Version: RStudio 2023.06.1+524 “Mountain Hydrangea”
Package Version:
Last updated: 23 February 2024
### --- Dataset importing and exporting -----
# Import in excel files
library(readxl)
# Export dataframes to excel files
library(writexl)
### ---- Data manipulation and visualization ----
## General data cleaning, manipulation, and visualization
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Multivariate data analysis
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
library(goeveg)
## This is GoeVeg 0.7.2 - build: 2024-02-06
## Multivariate visualization
library(ecodist)
##
## Attaching package: 'ecodist'
##
## The following object is masked from 'package:vegan':
##
## mantel
library(ecotraj)
## Loading required package: Rcpp
## Important functions to get the Bray-Curtis dissimilarity code to work
library(maditr)
##
## To modify variables or add new variables:
## let(mtcars, new_var = 42, new_var2 = new_var*hp) %>% head()
##
##
## Attaching package: 'maditr'
##
## The following objects are masked from 'package:dplyr':
##
## between, coalesce, first, last
##
## The following object is masked from 'package:purrr':
##
## transpose
##
## The following object is masked from 'package:readr':
##
## cols
## Making multipaneled figures
library(cowplot)
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(ggpubr)
##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:cowplot':
##
## get_legend
## Needed to draw images in plots
library(magick)
## Linking to ImageMagick 6.9.12.93
## Enabled features: cairo, fontconfig, freetype, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fftw, ghostscript, x11
## Prevents text from overlapping on figures
library(ggrepel)
## Ternary plots
library(ggtern)
## Registered S3 methods overwritten by 'ggtern':
## method from
## grid.draw.ggplot ggplot2
## plot.ggplot ggplot2
## print.ggplot ggplot2
## --
## Remember to cite, run citation(package = 'ggtern') for further info.
## --
##
## Attaching package: 'ggtern'
##
## The following objects are masked from 'package:ggplot2':
##
## aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
## ggsave, layer_data, theme_bw, theme_classic, theme_dark,
## theme_gray, theme_light, theme_linedraw, theme_minimal, theme_void
### ---- Model fitting and analysis ----
# Mixed-models
## Fits mixed models
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Functions to analyze significance and R^2 components of mixed-models
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
## Looks at variance components for mixed models
library(MuMIn)
## Anova function
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
# Fits negative binomial models
library(MASS)
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
# Post hoc analysis
library(emmeans)
## Below is code to install the pairwiseAdonis package from github
### library(devtools)
### install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")
library(pairwiseAdonis)
## Loading required package: cluster
# Indicator species analysis
library(indicspecies)
# Summary stats
##
library(mosaic)
## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2
##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.
##
## Attaching package: 'mosaic'
##
## The following objects are masked from 'package:car':
##
## deltaMethod, logit
##
## The following object is masked from 'package:lmerTest':
##
## rand
##
## The following object is masked from 'package:lme4':
##
## factorize
##
## The following object is masked from 'package:Matrix':
##
## mean
##
## The following object is masked from 'package:cowplot':
##
## theme_map
##
## The following objects are masked from 'package:goeveg':
##
## deg2rad, rad2deg
##
## The following object is masked from 'package:permute':
##
## shuffle
##
## The following objects are masked from 'package:dplyr':
##
## count, do, tally
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## stat
##
## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
## quantile, sd, t.test, var
##
## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sum
spring_cover <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Vegetative Cover Data 2019/Cleaned_Spring_Cover_Data_July_2019.xlsx")
fall_cover <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Vegetative Cover Data 2019/Cleaned_Fall_Cover_Data_September_2019.xlsx")
seed_rain <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Seed Rain 2019/Seed_Rain_Data_2019_FINAL.xlsx")
seed_bank <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Seed Bank 2020/Seed_Bank_2020_Cleaned.xlsx")
traits <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Traits and Species Info/SR_SB_VEG_Traits_Dataset.xlsx")
Study Sites
Our study sites were at Tucker Prairie Research Area (38○56’53.6” N, 91○59’40.0” W, Callaway County, MO) and at Prairie Fork Conservation Area (38○58’29.7” N, 91○44’03.3” W, Callaway County, MO). Tucker Prairie is a 59-hectare tract of unplowed North American tallgrass claypan prairie. Less than 0.5 % of intact tallgrass prairie ecosystem (i.e., never-been-plowed) remains in Missouri, and Tucker Prairie represents the last sizable remnant prairie in north central Missouri (Samson & Knopf, 1994). More than 250 species of plants inhabit Tucker Prairie, with representatives from 57 families and over 150 genera (Tropicos, n.d.). Native C4 grasses dominate the landscape including Sorghastrum nutans, Andropogon gerardii, Schizachyrium scoparium, and Sporobolus heterolepis (Drew, 1947; Rabinowitz & Rapp, 1980). Before 1957, Tucker Prairie was used for cattle grazing and occasional haying (Drew, 1947). From 1958 to 2002, Tucker Prairie was burned once every four years in the late winter or early spring (Rabinowitz & Rapp, 1980). Since 2002, Tucker Prairie has been managed on a 5-year burn rotation, where units are burned once in the late winter to early spring (Jan. – Mar.) and again 2-3 years later in the late summer to early fall (Aug. – Oct.).
Prairie Fork Conservation Area (PFCA) is over 450 hectares of former agricultural land being restored to tallgrass prairie and savanna ecosystems (Navarrete-Tindall et al., 2007; Newbold et al., 2019). From 2004 onwards, 16-25 hectares are newly seeded each year with native prairie species collected from Tucker Prairie and other nearby remnant prairies (Newbold et al., 2019). As a result, the plots within PFCA represent a chronosequence of reconstructed tallgrass prairies that are mainly comprised of Tucker Prairie descendants. Prior to seeding, glyphosate-resistant crops (i.e., corn and soybeans) are planted and harvested for at least three years to reduce the existing weed seed bank (Newbold et al., 2019). A seed mix containing an average of 179 species of native forbs, legumes, grasses, rushes, and sedges are broadcast seeded at a rate of 13.4 to 18.2 kg/ha once during the dormant season (Newbold et al., 2019). All mixes are comprised of at least 75% of the same species each year that are representative of nearby remnant prairies such as Tucker Prairie (Newbold et al., 2019). Reconstructions are managed using a 2-4-year burn schedule and annual spraying of invasive species (Newbold et al., 2019). For additional details on PFCA management, see Newbold et al. 2019. To capture changes in seed rain dynamics during the restoration process, we grouped our reconstructed sites into three categories, as defined by Newbold et al. (2019) as being compositionally different. We measured the seed rain in an old reconstruction (seeded in 2004), middle aged reconstruction (seeded in 2012 and 2013), and a young reconstruction (seeded in 2017).
Seed rain sampling
In May 2019, we deployed artificial turf grass seed traps (0.1 x 0.1 m) in Tucker Prairie and at each of the focal PFCA restorations. We used artificial turf grass traps instead of sticky traps due to their durability and resistance to freezing (Molau & Per Mølgaard, 1996; Bass, 2015). Turfgrass traps also do not lose effectiveness over time, like sticky traps, which are prone to contamination by soil, non-seed plant debris, and dead insects, making seed recovery difficult (Arruda et al., 2020). Furthermore, unlike funnel traps, they are easy to install and collect with minimal disturbance to surrounding vegetation and soil and do not fill with water (Schott, 1995). At each site, we randomly established ten, 5 m long transects, placed a minimum of 50 m apart from each other. We then placed traps at 1 m intervals along each transect and affixed them to the ground using ground staples (n = 5 traps per transect, 50 traps per site). We collected and replaced seed traps every 2-weeks during the growing season from May 31st to December 12th, 2019. Hereafter we refer to all captured diaspores as “seeds” despite catching a variety of fruit types (e.g., achene). After collection, all traps from a transect and sampling period were grouped together and sieved through a series of soil sieves (1 mm, 500 mm, 250 mm mesh). From these layers, we then picked out and identified seeds to the lowest taxonomic level possible using a reference collection of species inhabiting our study areas and identification guides (Coons et al., n.d.; A. C. Martin & Barkley, 1961).
Because of their extremely small size (0.3 – 0.5 mm; Yatskievych, 1999), we estimated the number of rush seeds (Juncus sp.) when there were over 200 seeds present in a sample. We first sieved the samples through 180 mm and 150 mm mesh soil sieves. Then we subsampled each sieved layer by calculating the average number of rush seeds per 1 cm2 area (n = 3) and multiplying that average by the total area covered by the sample layer. Lastly, we summed the number of estimated rush seeds per layer to calculate the total number of rush seeds in a sample.
Vegetation sampling
In 2019, at the same transects, we sampled species cover twice, once in the early summer (June - July) and again at peak biomass (September). During the first sampling period, we only recorded species presence. At peak biomass, we also sampled percent aerial species cover in a 1m2 area centered around each seed rain trap. We identified species according to Yatskievych (1999).
Seed bank sampling
In March 2020, before the growing season, we collected soil samples to assess the soil seed bank at each focal prairie. To compare the soil seed bank and the previous year’s seed rain, we collected 5 (10 x 10 x 10 cm3; 1000 cm3) soil cores approximately 0.5 m away from where we captured seed rain at each transect in 2019. We allowed the soil samples to dry before removing all non-seed plant material (e.g., roots and rhizomes) to ensure seedlings only germinated from seeds. We then homogenized all samples collected from a particular transect and subsampled ~1500 cm3 of soil to spread ~1 cm deep over sterile potting soil in plastic germination trays. Starting July 2020, we placed the trays (n = 40) in a greenhouse and watered them when dry. Amongst the 40 trays containing soil samples, we also placed an additional 10 trays containing only sterile potting soil to assess whether any contamination occurred from external sources. Once identifiable, we removed seedlings from trays. After germination ceased in July 2021, we placed all trays into a cold room for ~4 months to replicate conditions necessary to break dormancy for any still dormant seeds. Post-vernalization, we returned the trays to the greenhouse, stirred the soil samples, and monitored for any additional germination. We ended the study one year post-vernalization.
# Make the six letter species codes for the spring and fall cover uppercase
spring_cover[[3]] <- toupper(spring_cover[[3]])
fall_cover[[4]] <- toupper(fall_cover[[4]])
# Remove unnecessary columns from datasets
### Spring cover - Only keep Site, Transect, SPP6, Scientific Name
spring_cover <- spring_cover[,-c(5,6)]
### Fall cover - Keep everything
### Seed rain - Only keep Plot, Transect, Sampling Date, Week, SPP6, Number, Unknown
seed_rain <- seed_rain[, -c(7,8)]
### Seed bank - Keep everything
# Change names of columns to better reflect datasets (e.g. seed rain: number -> Number_Seeds)
colnames(seed_bank) <- c("Site", "Transect", "SPP6", "Number_Seedlings")
colnames(seed_rain) <- c("Site", "Transect", "Sampling_Date", "Week", "SPP6", "Number_Seeds", "Unknown")
# Change transect PFCA 2 - 10A to just 10 for the cover datasets
spring_cover <- spring_cover %>%
mutate_if(is.character,
str_replace_all, pattern = "10A", replacement = "10")
fall_cover <- fall_cover %>%
mutate_if(is.character,
str_replace_all, pattern = "10a", replacement = "10")
# Make Transect and Plot variables factors for all datasets
## Spring cover
spring_cover$Transect <- as.factor(spring_cover$Transect)
## Fall cover
fall_cover$Transect <- as.factor(fall_cover$Transect)
fall_cover$Plot <- as.factor(fall_cover$Plot)
## Seed Bank
seed_bank$Transect <- as.factor(seed_bank$Transect)
## Seed Rain
seed_rain$Transect <- as.factor(seed_rain$Transect)
####### Clean Spring
# 0.) Make a frequency column
spring_cover$Cover <- rep(1, nrow(spring_cover))
####### Clean Fall
# 0.) Make anything less than 1 = 1
fall_cover <- fall_cover %>%
mutate(Cover = ifelse(Cover < 1, 1, Cover))
# 1.) get a sum per species per transect of cover.
fall_sum <- fall_cover %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarize(Cover = sum(Cover))
## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.
fall_sum
## # A tibble: 1,215 × 5
## # Groups: Site, Transect, SPP6 [1,214]
## Site Transect SPP6 Scientific_Name Cover
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ACAVIR Acalypha virginica 11.5
## 2 PFCA 1 1 AMBART Ambrosia artemisiifolia 6
## 3 PFCA 1 1 BARVUL Barbarea vulgaris 3.5
## 4 PFCA 1 1 CHAFAS Chamaecrista fasciculata 388
## 5 PFCA 1 1 CONCAN Conyza canadensis 1
## 6 PFCA 1 1 ERISPP Erigeron spp. 2
## 7 PFCA 1 1 ERIVIL Eriochloa villosa 2
## 8 PFCA 1 1 EUPSER Eupatorium serotinum 3
## 9 PFCA 1 1 KUMSTI Kummerowia stipulacea 6
## 10 PFCA 1 1 LESCAP Lespedeza capitata 1
## # ℹ 1,205 more rows
####### Seed Rain
# 1.) get a sum per species per transect of cover.
seed_rain_sum <- seed_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
# Clean in the same way you did the seed rain analysis
for(i in 1:nrow(seed_rain_sum)){
if(seed_rain_sum[i,3] == "AGAFAS"){seed_rain_sum[i,3] <- "AGASPP"}
if(seed_rain_sum[i,3] == "EUPSER"){seed_rain_sum[i,3] <- "EUPSPP"}
if(seed_rain_sum[i,3] == "GENPUB"){seed_rain_sum[i,3] <- "GENSPP"}
}
####### Seed Bank
# 1.) get a sum per species per transect of cover.
seed_bank_sum <- seed_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
####### Remove Spring Unknowns
spring_cover <- filter(spring_cover, !grepl('UNK', SPP6))
# Looks good
list(unique(spring_cover$SPP6))
## [[1]]
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART" "AMBPSI"
## [9] "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS" "BAPBRA"
## [17] "BARVUL" "BIDARI" "BLECIL" "BROJAP" "CAMRAD" "CAPBUR" "CARSPP" "CARBIC"
## [25] "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL" "CORSPP"
## [33] "CORTRI" "CROSAG" "CYPACU" "CYPECH" "DALPUR" "DESILL" "DESSES" "DESSPP"
## [41] "DICLAN" "DIGISC" "ECHPAL" "ELETEN" "ELYVIR" "ERISPP" "ERYYUC" "EUPCOR"
## [49] "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENPUB"
## [57] "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HYPPUN" "JUNBRA" "JUNMAR"
## [65] "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPVIR" "LESCAP" "LESCUN"
## [73] "LESVIR" "LIAPYC" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
## [81] "MELSPP" "MONFIS" "MUHSPP" "MYOVER" "OENBIE" "OXADIL" "PARINT" "PENDIG"
## [89] "PLAVIR" "POAPRA" "POLSAN" "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO"
## [97] "RATPIN" "RHUCOP" "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI"
## [105] "SALAZU" "SCHSCO" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [113] "SISSPP" "SOLALT" "SOLCAR" "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOHET"
## [121] "STRLEI" "SYMERI" "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB"
## [129] "SYMPRA" "SYMSPP" "THLARV" "OENSPP" "TRAOHI" "TRIPER" "TRIREP" "TRISPP"
## [137] "ULMSPP" "SCISPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [145] "ZIZAUR" "SCLTRI"
####### Remove Fall Unknowns
fall_sum <- filter(fall_sum, !grepl('UNK', SPP6))
# Looks good
list(unique(fall_sum$SPP6))
## [[1]]
## [1] "ACAVIR" "AMBART" "BARVUL" "CHAFAS" "CONCAN" "ERISPP" "ERIVIL" "EUPSER"
## [9] "KUMSTI" "LESCAP" "PENDIG" "PLAVIR" "RUDHIR" "SETFAB" "SOLALT" "STRLEI"
## [17] "SYMPIL" "ACHMIL" "CORTRI" "DALPUR" "DIGISC" "EUTGYM" "JUNSPP" "PYCTEN"
## [25] "RATPIN" "RUDSUB" "SCISPP" "SOLSPE" "SORNUT" "SYMERI" "SYMNOV" "SYMORB"
## [33] "TRICAM" "BROJAP" "HELMOL" "HYPPUN" "OENBIE" "SOLRIG" "SYMLAN" "AGRHYE"
## [41] "KUMSTR" "LESCUN" "LIAPYC" "PARINT" "SOLCAR" "TAROFF" "BAPALB" "ERASPE"
## [49] "JUNVIR" "LESVIR" "SYMPRA" "VERARV" "CERSPP" "ECHCRU" "OXADIL" "PLAMAJ"
## [57] "SILINT" "SILLAC" "VERSPP" "CORLAN" "CYPECH" "CYPSTR" "EUPCOR" "MEDLUP"
## [65] "SCHSCO" "BIDARI" "FESPAR" "MONFIS" "RANABO" "SOLNEM" "OENFIL" "ULMPUM"
## [73] "DESSPP" "ELYCAN" "AGATEN" "BAPBRA" "CARSPP" "CIRALT" "ERYYUC" "GENAND"
## [81] "KOEMAC" "LACSER" "PYCPIL" "RUMCRI" "SPOHET" "SYMOBL" "VERHEL" "BLECIL"
## [89] "BOLAST" "CYPACU" "SYMLAT" "SYMOOL" "ECHPAL" "MELSPP" "SYMANO" "SYMLAE"
## [97] "CORPAL" "ANDGER" "EUPALT" "POAPRA" "ALOCAR" "LIASPP" "LYTALA" "SPOCOM"
## [105] "SPILAC" "PANVIR" "SENMAR" "TRIFLA" "HELSPP" "ELYVIR" "HELHEL" "MYOVER"
## [113] "SALAZU" "TRAOHI" "BAPAUS" "LESPRO" "DESSES" "ROSCAR" "AMOCAN" "LEPDEN"
## [121] "CAMRAD" "RUBSPP" "ACERUB" "ZIZAUR" "DICLAN" "EUPPER" "GALOBT" "LYSLAN"
## [129] "POTSIM" "SETPAR" "SOLMIS" "VIOSAG" "AGAFAS" "ANTNEG" "CROSAG" "LINSUL"
## [137] "RHUGLA" "FRAVIR" "GENPUB" "MUHGLA" "JUNBRA" "ULMSPP" "POLVER" "TRISPP"
## [145] "CORSPP" "POLSAN" "SILANT" "TRIPER" "RHUCOP"
####### Remove Seed Rain Unknowns
# 1.) get a sum per species per transect of cover.
seed_rain_sum <- filter(seed_rain_sum, !grepl('UNK', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('PANICUM?', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('NONE', SPP6))
# Looks good
list(unique(seed_rain_sum$SPP6))
## [[1]]
## [1] "ACAVIR" "AMBART" "ANAMIN" "BOLAST" "CERSPP" "CHAFAS" "CHEALB"
## [8] "CONCAN" "DIGISC" "ERISPP" "ERIVIL" "EUPSPP" "KUMSPP" "LEPVIR"
## [15] "MEDLUP" "OXADIL" "PENDIG" "PLAVIR" "PYCTEN" "RUDHIR" "SETSPP"
## [22] "SOLALT" "SORNUT" "SYMSPP" "THLARV" "TRIPER" "VEROSP" "ALOCAR"
## [29] "CORTRI" "CYPECH" "DICLAN" "ERYYUC" "HYPPUN" "MONFIS" "MYOVER"
## [36] "OENBIE" "RATPIN" "SCIGEO" "SOLRIG" "SPHOBT" "STRLEI" "TRIFLA"
## [43] "ACHMIL" "BARVUL" "FESARU" "GERCAR" "LESCUN" "LIAPYC" "PANCAP"
## [50] "SYMNOV" "BIDARI" "CARBUS" "DESSPP" "ERASPE" "GALAPA" "LESVIR"
## [57] "POAPRA" "RUDSUB" "VERHAS" "VERSPP" "ANDGER" "CARFES" "ECHCRU"
## [64] "EREHIE" "PLAOCC" "SILINT" "AGASPP" "BROJAP" "CAPBUR" "CIRALT"
## [71] "CORLAN" "MELSPP" "MOLVER" "SCHSCO" "SCLTRI" "CYPSTR" "FESPAR"
## [78] "JUNSPP" "PYCPIL" "TAROFF" "AMATUB" "LESCAP" "OENFIL" "KOEMAC"
## [85] "DAUCAR" "LEUVUL" "VULOCT" "CYPACU" "HORPUS" "SCIPEN" "ECHPAL"
## [92] "SOLNEM" "CORPAL" "LINSUL" "DIAARM" "HELMOL" "PERLON" "RUMCRI"
## [99] "BLECIL" "LYTALA" "SPOHET" "ELESPP1" "SPOCOM" "CARCEP" "AGRHYE"
## [106] "GENSPP" "HELHEL" "TRAOHI" "BAPALB" "CARANN" "IPOHED" "AMOCAN"
## [113] "SALAZU" "CARBIC" "EUPCOR" "EUTGYM" "GALOBT" "RUBSPP" "SETPAR"
## [120] "ELESPP" "LYSLAN" "POLSPP" "HYPHIR" "LYCAME" "SILANT" "CARDSP"
## [127] "CROSAG" "LOBSPI" "VIOSAG"
####### Remove Seed Bank Unknowns
seed_bank_sum <- filter(seed_bank_sum, !grepl('UNK', SPP6))
# Looks good
list(sort(unique(seed_bank_sum$SPP6)))
## [[1]]
## [1] "ABUTHE" "ACAVIR" "ACHMIL" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "ANDGER"
## [9] "BARVUL" "CARFES" "CARHIR" "CARPAR" "CERGLO" "CHAFAS" "CIRALT" "CONCAN"
## [17] "CORLAN" "CROSAG" "CYPACU" "CYPSPP" "CYPSTR" "DESSPP" "DICLAN" "DIGISC"
## [25] "DIGSAN" "DIPLAC" "ERASPE" "ERIANN" "ERISTR" "EUPHUM" "EUPSER" "EUTGYM"
## [33] "GERCAR" "HORPUS" "HYPPUN" "IPOHED" "JUNINT" "JUNMAR" "KUMSTI" "KUMSTR"
## [41] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LYSLAN" "MEDLUP" "MELSPP" "MOLVER"
## [49] "OENSPP" "OXADIL" "PANCAP" "PANDIC" "PENDIG" "PERLON" "PLAOCC" "PLAPUS"
## [57] "PLAVIR" "POACHA" "POAPRA" "POPDEL" "POROLE" "PYCPIL" "PYCTEN" "RATPIN"
## [65] "RUDHIR" "RUMCRI" "SCHSCO" "SETFAB" "SETGLA" "SILANT" "SOLALT" "SOLCAR"
## [73] "SORNUT" "SPHOBT" "STRLEI" "SYMLAE" "SYMPIL" "TAROFF" "THLARV" "TOXRAD"
## [81] "TRAOHI" "TRIFLA" "TRIPER" "TRIREP" "VERHAS" "VERPER" "VERTHA"
# Full join the two cover data sets
cover_merged <- full_join(spring_cover, fall_sum, by = c("Site", "Transect", "SPP6", "Scientific_Name", "Cover")) %>%
mutate_all(~replace(., is.na(.), 0))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
# Make sure there is no weirdness in the join
#View(cover_merged)
### Group by and then select the SPP6 with the max cover between the two datasets
### Since Spring only had presence/absence we don't want to inflate the fall cover artificially if we saw that species again
cover_merged_max <- cover_merged %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarise(Cover=max(Cover))
## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.
## Not exactly sure what to do about SYMSPP and KUMSPP (things we couldn't ID in the spring but could in the fall)
### Going to drop SYMSPP from the dataset since every plot that had SYMSPP in the spring had identified Symphyotrichum in the fall
cover_merged_max <- filter(cover_merged_max , !grepl('SYMSPP', SPP6))
### Checked KUMSPP only PFCA 2 - 3 and PFCA 2 - 9 did not see Kummerowia again in the fall; will drop KUMSPP for all other plots
### Removes all the PFCA 1 KUMSPPs
cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 1"),]
### Removes all PFCA 2 KUMSPPs that were not in PFCA 2 - 3 and PFCA 2- 9
cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 2" & cover_merged_max$Transect%in% c(2,4,5,8)),]
### Check all the species codes
sort(unique(cover_merged_max$SPP6))
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGAFAS" "AGATEN" "AGRGIG" "AGRHYE" "ALOCAR"
## [9] "AMBART" "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
## [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
## [25] "CAPBUR" "CARBIC" "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN"
## [33] "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR"
## [41] "DALPUR" "DESILL" "DESSES" "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
## [49] "ELETEN" "ELYCAN" "ELYVIR" "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPALT"
## [57] "EUPCOR" "EUPPER" "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA"
## [65] "GALOBT" "GENAND" "GENPUB" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP"
## [73] "HYPPUN" "JUNBRA" "JUNMAR" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "KUMSTI"
## [81] "KUMSTR" "LACSER" "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR"
## [89] "LIAPYC" "LIASPP" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
## [97] "MELSPP" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP"
## [105] "OXADIL" "PANVIR" "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSAN"
## [113] "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMANO" "SYMERI" "SYMLAE"
## [153] "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB" "SYMPIL" "SYMPRA"
## [161] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [169] "ULMPUM" "ULMSPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [177] "ZIZAUR"
### Cover grouping so that it works with the seed rain dataset
# Code to lump certain species together
cover_to_merge_rain <- cover_merged_max
for(i in 1:nrow(cover_to_merge_rain)){
# Group all Symphyotrichum except SYMNOV
## SYMPIL
## SYMLAN
## SYMLAT
## SYMERI
## SYMANO
## SYMLAE
## SYMPRA
## SYMOBL
## SYMOOL
if(cover_to_merge_rain[i,3] == "SYMPIL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAN"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAT"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMERI"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMANO"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAE"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMPRA"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMOBL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMOOL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
# Group all Carex
## CARBIC
if(cover_to_merge_rain[i,3] == "CARBIC"){cover_to_merge_rain[i,3] <- "CARSPP"}
# Group all Kummerowia
## KUMSTI
## KUMSTR
if(cover_to_merge_rain[i,3] == "KUMSTI"){cover_to_merge_rain[i,3] <- "KUMSPP"}
if(cover_to_merge_rain[i,3] == "KUMSTR"){cover_to_merge_rain[i,3] <- "KUMSPP"}
# Group all Juncus
## JUNMAR
## JUNBBRA
if(cover_to_merge_rain[i,3] == "JUNMAR"){cover_to_merge_rain[i,3] <- "JUNSPP"}
if(cover_to_merge_rain[i,3] == "JUNBRA"){cover_to_merge_rain[i,3] <- "JUNSPP"}
# Group all Agalinis
## AGATEN
## AGAFAS
if(cover_to_merge_rain[i,3] == "AGATEN"){cover_to_merge_rain[i,3] <- "AGASPP"}
if(cover_to_merge_rain[i,3] == "AGAFAS"){cover_to_merge_rain[i,3] <- "AGASPP"}
# Group all Eleocharis
## ELETEN
if(cover_to_merge_rain[i,3] == "ELETEN"){cover_to_merge_rain[i,3] <- "ELESPP"}
# Group all Gentian
## GENPUB
## GENAND
if(cover_to_merge_rain[i,3] == "GENPUB"){cover_to_merge_rain[i,3] <- "GENSPP"}
if(cover_to_merge_rain[i,3] == "GENAND"){cover_to_merge_rain[i,3] <- "GENSPP"}
# Group all Setaria except Setaria parviflora
## SETGLA
## SETFAB
if(cover_to_merge_rain[i,3] == "SETGLA"){cover_to_merge_rain[i,3] <- "SETSPP"}
if(cover_to_merge_rain[i,3] == "SETFAB"){cover_to_merge_rain[i,3] <- "SETSPP"}
# Group all Veronica
## VERARV
if(cover_to_merge_rain[i,3] == "VERARV"){cover_to_merge_rain[i,3] <- "VEROSP"}
# Group all Desmodium
## DESSES
if(cover_to_merge_rain[i,3] == "DESSES"){cover_to_merge_rain[i,3] <- "DESSPP"}
# Group all Eupatorium
## EUPALT
## EUPSER
## EUPPER
if(cover_to_merge_rain[i,3] == "EUPALT"){cover_to_merge_rain[i,3] <- "EUPSPP"}
if(cover_to_merge_rain[i,3] == "EUPSER"){cover_to_merge_rain[i,3] <- "EUPSPP"}
if(cover_to_merge_rain[i,3] == "EUPPER"){cover_to_merge_rain[i,3] <- "EUPSPP"}
# Group all Polygala
## POLVER
## POLSAN
if(cover_to_merge_rain[i,3] == "POLVER"){cover_to_merge_rain[i,3] <- "POLSPP"}
if(cover_to_merge_rain[i,3] == "POLSAN"){cover_to_merge_rain[i,3] <- "POLSPP"}
}
### Check that it changed the codes
sort(unique(cover_to_merge_rain$SPP6))
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART"
## [9] "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS"
## [17] "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD" "CAPBUR"
## [25] "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL"
## [33] "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DALPUR" "DESILL"
## [41] "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ELYCAN" "ELYVIR"
## [49] "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
## [57] "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELFLE" "HELHEL"
## [65] "HELMOL" "HELSPP" "HYPPUN" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER"
## [73] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LIAPYC" "LIASPP"
## [81] "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP" "MONFIS"
## [89] "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP" "OXADIL" "PANVIR"
## [97] "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSPP" "POTSIM" "PYCPIL"
## [105] "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA" "ROSCAR" "ROSSPP" "RUBSPP"
## [113] "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO" "SCISPP" "SCLTRI" "SENMAR"
## [121] "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC" "SISSPP" "SOLALT" "SOLCAR"
## [129] "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT" "SPHOBT" "SPILAC" "SPOCOM"
## [137] "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP" "TAROFF" "THLARV" "TRAOHI"
## [145] "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP" "ULMPUM" "ULMSPP" "VERHEL"
## [153] "VEROSP" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT" "ZIZAUR"
### Sum everything by species by transect again
cover_to_merge_rain <- cover_to_merge_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Cover = sum(Cover))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
### Seed grouping so that it works with the cover dataset
seed_rain_names <- left_join(seed_rain_sum, traits)
## Joining with `by = join_by(SPP6)`
rain_to_merge_cover <- seed_rain_names[,c(1,2,3,4,5)]
for(i in 1:nrow(rain_to_merge_cover)){
# Group all Carex
## CARBUS
## CARBIC
## CARANN
## CARFES
## CARCEP
if(rain_to_merge_cover[i,3] == "CARBUS"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARBIC"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARANN"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARFES"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARCEP"){rain_to_merge_cover[i,3] <- "CARSPP"}
# Group all Eleocharis
## ELESPP1
if(rain_to_merge_cover[i,3] == "ELESPP1"){rain_to_merge_cover[i,3] <- "ELESPP"}
# Group all Scirpus
if(rain_to_merge_cover[i,3] == "SCIGEO"){rain_to_merge_cover[i,3] <- "SCISPP"}
if(rain_to_merge_cover[i,3] == "SCIPEN"){rain_to_merge_cover[i,3] <- "SCISPP"}
}
### Sum everything by species by transect again
rain_to_merge_cover <- rain_to_merge_cover %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
# Looks good
sort(unique(rain_to_merge_cover$SPP6))
## [1] "ACAVIR" "ACHMIL" "AGASPP" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "AMOCAN"
## [9] "ANAMIN" "ANDGER" "BAPALB" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP"
## [17] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "CONCAN"
## [25] "CORLAN" "CORPAL" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DAUCAR"
## [33] "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ERASPE"
## [41] "EREHIE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
## [49] "FESPAR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELHEL" "HELMOL" "HORPUS"
## [57] "HYPHIR" "HYPPUN" "IPOHED" "JUNSPP" "KOEMAC" "KUMSPP" "LEPVIR" "LESCAP"
## [65] "LESCUN" "LESVIR" "LEUVUL" "LIAPYC" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN"
## [73] "LYTALA" "MEDLUP" "MELSPP" "MOLVER" "MONFIS" "MYOVER" "OENBIE" "OENFIL"
## [81] "OXADIL" "PANCAP" "PENDIG" "PERLON" "PLAOCC" "PLAVIR" "POAPRA" "POLSPP"
## [89] "PYCPIL" "PYCTEN" "RATPIN" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU"
## [97] "SCHSCO" "SCISPP" "SCLTRI" "SETPAR" "SETSPP" "SILANT" "SILINT" "SOLALT"
## [105] "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV"
## [113] "SYMSPP" "TAROFF" "THLARV" "TRAOHI" "TRIFLA" "TRIPER" "VERHAS" "VEROSP"
## [121] "VERSPP" "VIOSAG" "VULOCT"
# Drop the scientific names they are causing trouble with the join
cover_to_merge_rain <- cover_to_merge_rain
cover_rain <- full_join(cover_to_merge_rain, rain_to_merge_cover) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:
# Drop TP transect 7 because of missing data
cover_rain <- cover_rain[!(cover_rain$Site%in%"TP" & cover_rain$Transect%in% "7"),]
### Cover grouping so that it works with the seed bank dataset
# Code to lump certain species together
cover_to_merge_bank <- cover_merged_max
for(i in 1:nrow(cover_to_merge_bank)){
# Group all Desmodium together
## DESSES
if(cover_to_merge_bank[i,3] == "DESSES"){cover_to_merge_bank[i,3] <- "DESSPP"}
# Group all Cyperus together
## CYPSTR
## CYPECH
## CYPACU
if(cover_to_merge_bank[i,3] == "CYPSTR"){cover_to_merge_bank[i,3] <- "CYPSPP"}
if(cover_to_merge_bank[i,3] == "CYPECH"){cover_to_merge_bank[i,3] <- "CYPSPP"}
if(cover_to_merge_bank[i,3] == "CYPACU"){cover_to_merge_bank[i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(cover_to_merge_bank[i,3] == "OENFIL"){cover_to_merge_bank[i,3] <- "OENSPP"}
if(cover_to_merge_bank[i,3] == "OENBIE"){cover_to_merge_bank[i,3] <- "OENSPP"}
}
cover_to_merge_bank <- cover_to_merge_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Cover = sum(Cover))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
### Seed Bank grouping so that it works with the cover dataset
bank_to_merge_cover <- seed_bank_sum
# Code to lump certain species together
for(i in 1:nrow(bank_to_merge_cover)){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_cover[i,3] == "CERGLO"){bank_to_merge_cover[i,3] <- "CERSPP"}
# Group all Cyperus together
## CYPSTR
## CYPACU
if(bank_to_merge_cover[i,3] == "CYPSTR"){bank_to_merge_cover[i,3] <- "CYPSPP"}
if(bank_to_merge_cover[i,3] == "CYPACU"){bank_to_merge_cover[i,3] <- "CYPSPP"}
# Group all the Carex together
if(bank_to_merge_cover[i,3] == "CARFES"){bank_to_merge_cover[i,3] <- "CARSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_cover[i,3] == "ERIANN"){bank_to_merge_cover[i,3] <- "ERISPP"}
if(bank_to_merge_cover[i,3] == "ERISTR"){bank_to_merge_cover[i,3] <- "ERISPP"}
# Change Juncus interior to Juncus sp. (since we were able to ID JUNMAR in the cover but not JUNINT)
if(bank_to_merge_cover[i,3] == "JUNINT"){bank_to_merge_cover[i,3] <- "JUNSPP"}
}
bank_to_merge_cover <- bank_to_merge_cover %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
# Drop the scientific names they are causing trouble with the join
cover_bank <- full_join(cover_to_merge_bank, bank_to_merge_cover) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:
# Drop TP transect 7 because of missing fall cover data
cover_bank <- cover_bank[!(cover_bank$Site%in%"TP" & cover_bank$Transect%in% "7"),]
### Seed Bank grouping so that it works with the cover dataset
bank_to_merge_rain <- seed_bank_sum
# Code to lump certain species together
for(i in 1:nrow(bank_to_merge_rain)){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_rain[i,3] == "CERGLO"){bank_to_merge_rain[i,3] <- "CERSPP"}
# Group all Cyperus together
## CYPSTR
## CYPACU
if(bank_to_merge_rain[i,3] == "CYPSTR"){bank_to_merge_rain[i,3] <- "CYPSPP"}
if(bank_to_merge_rain[i,3] == "CYPACU"){bank_to_merge_rain[i,3] <- "CYPSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_rain[i,3] == "ERIANN"){bank_to_merge_rain[i,3] <- "ERISPP"}
if(bank_to_merge_rain[i,3] == "ERISTR"){bank_to_merge_rain[i,3] <- "ERISPP"}
# Group all the Juncus together
## JUNINT
## JUNMAR
if(bank_to_merge_rain[i,3] == "JUNINT"){bank_to_merge_rain[i,3] <- "JUNSPP"}
if(bank_to_merge_rain[i,3] == "JUNMAR"){bank_to_merge_rain[i,3] <- "JUNSPP"}
# Group all the Setaria together
if(bank_to_merge_rain[i,3] == "SETFAB"){bank_to_merge_rain[i,3] <- "SETSPP"}
if(bank_to_merge_rain[i,3] == "SETGLA"){bank_to_merge_rain[i,3] <- "SETSPP"}
# Group all the Symphyotrichum together
## SYMPIL
## SYMLAE
if(bank_to_merge_rain[i,3] == "SYMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
if(bank_to_merge_rain[i,3] == "SYMLAE"){bank_to_merge_rain[i,3] <- "SYMSPP"}
if(bank_to_merge_rain[i,3] == "SIMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
# Group all the Kummerowia together
## KUMSTR
## KUMSTI
if(bank_to_merge_rain[i,3] == "KUMSTR"){bank_to_merge_rain[i,3] <- "KUMSPP"}
if(bank_to_merge_rain[i,3] == "KUMSTI"){bank_to_merge_rain[i,3] <- "KUMSPP"}
# Group all the Veronica together
## VERPER
if(bank_to_merge_rain[i,3] == "VERPER"){bank_to_merge_rain[i,3] <- "VEROSP"}
# Group all the Eupatorium together
## EUPSER
if(bank_to_merge_rain[i,3] == "EUPSER"){bank_to_merge_rain[i,3] <- "EUPSPP"}
# Merge all the Cardamine together
## CARHIR
## CARPAR
if(bank_to_merge_rain[i,3] == "CARHIR"){bank_to_merge_rain[i,3] <- "CARDSP"}
if(bank_to_merge_rain[i,3] == "CARPAR"){bank_to_merge_rain[i,3] <- "CARDSP"}
}
bank_to_merge_rain <- bank_to_merge_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
rain_to_merge_bank <- seed_rain_sum
for(i in 1:nrow(rain_to_merge_bank)){
# Group all Cyperus together
## CYPSTR
## CYPECH
## CYPACU
if(rain_to_merge_bank[i,3] == "CYPSTR"){rain_to_merge_bank[i,3] <- "CYPSPP"}
if(rain_to_merge_bank[i,3] == "CYPECH"){rain_to_merge_bank[i,3] <- "CYPSPP"}
if(rain_to_merge_bank[i,3] == "CYPACU"){rain_to_merge_bank[i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(rain_to_merge_bank[i,3] == "OENFIL"){rain_to_merge_bank[i,3] <- "OENSPP"}
if(rain_to_merge_bank[i,3] == "OENBIE"){rain_to_merge_bank[i,3] <- "OENSPP"}
}
rain_to_merge_bank <- rain_to_merge_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))
## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.
### Merge the seed rain and seed bank datasets together
rain_bank <- full_join(bank_to_merge_rain, rain_to_merge_bank) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:
Wasn’t in the cover at all but we caught seeds
Didn’t catch seeds but was in the cover
Notable (> 10% cover observed)
Minor ( 1% < cover < 10% observed)
Trace (cover 1% observed)
## # A tibble: 10,959 × 5
## # Groups: Site, Transect [39]
## Site Transect SPP6 Abundance Type
## <chr> <fct> <chr> <dbl> <chr>
## 1 PFCA 1 1 ACAVIR 11.5 Aboveground
## 2 PFCA 1 10 ACAVIR 4 Aboveground
## 3 PFCA 1 2 ACAVIR 4.5 Aboveground
## 4 PFCA 1 3 ACAVIR 2 Aboveground
## 5 PFCA 1 4 ACAVIR 4 Aboveground
## 6 PFCA 1 5 ACAVIR 3 Aboveground
## 7 PFCA 1 6 ACAVIR 2.5 Aboveground
## 8 PFCA 1 7 ACAVIR 3 Aboveground
## 9 PFCA 1 8 ACAVIR 1 Aboveground
## 10 PFCA 1 9 ACAVIR 5 Aboveground
## # ℹ 10,949 more rows
# Code obtained and modified from Lauren and Anu's seed rain vs. vegetation paper
##### run dissimilarity loops
Cover_rain_dis_results <-matrix(nrow=0,ncol=7)
sites <-as.vector(unique(cover_rain_dis$Site))
for(s in 1:length(sites)){
temp <-subset(cover_rain_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
#for(b in 1:length(blocks)){
#temp_b<-subset(temp, block==blocks[b])
#plots<-as.vector(unique(temp_b$plot))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"hellinger"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
trans_int <- cbind(plot_cast$Type,trunc(plot_cast[,-1])) #integers only
names(trans_int)[1]<-"trt"
trans_int_norm <- cbind(plot_cast$Type,trunc(decostand(plot_cast[,-1], "normalize")*100)) #integers only and normalize
names(trans_int)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
distance_kulczynski<-vegdist(trans_tot[,-1], method = "kulczynski")
distance_cao<-vegdist(trans_int_norm[,-1], method = "cao")
distance_morisita<-vegdist(trans_int[,-1], method = "morisita")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1],
distance_cao[1], distance_morisita[1], distance_kulczynski[1])
Cover_rain_dis_results <-rbind(Cover_rain_dis_results, new.row)
}
print(s/length(sites))
}
## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1
row.names(Cover_rain_dis_results)<-NULL
Cover_rain_dis_results<- data.frame(Cover_rain_dis_results)
names(Cover_rain_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard",
"dissimilarity_cao", "dissimilarity_morisita", "dissimilarity_kulczynski")
##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = Cover_rain_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.148341 -0.033152 0.005915 0.046429 0.104723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.51363 0.01974 26.018 < 2e-16 ***
## SitePFCA 1 -0.13402 0.02721 -4.925 2.01e-05 ***
## SitePFCA 2 -0.07069 0.02721 -2.598 0.0136 *
## SitePFCA 3 0.02310 0.02721 0.849 0.4017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05922 on 35 degrees of freedom
## Multiple R-squared: 0.5505, Adjusted R-squared: 0.512
## F-statistic: 14.29 on 3 and 35 DF, p-value: 3.056e-06
## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 2.37434 1 676.96 < 2.2e-16 ***
## Site 0.15036 3 14.29 3.056e-06 ***
## Residuals 0.12276 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cover_rain.mod.bray)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.13402090 -0.207406554 -0.060635237 0.0001142
## PFCA 2-TP -0.07069333 -0.144078991 0.002692326 0.0624794
## PFCA 3-TP 0.02310178 -0.050283878 0.096487439 0.8306569
## PFCA 2-PFCA 1 0.06332756 -0.008100791 0.134755919 0.0975147
## PFCA 3-PFCA 1 0.15712268 0.085694321 0.228551031 0.0000055
## PFCA 3-PFCA 2 0.09379511 0.022366758 0.165223468 0.0060381
##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = Cover_rain_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.104048 -0.037405 -0.006949 0.037294 0.192794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.54155 0.02038 26.578 < 2e-16 ***
## SitePFCA 1 -0.06193 0.02809 -2.205 0.03412 *
## SitePFCA 2 -0.07926 0.02809 -2.822 0.00782 **
## SitePFCA 3 0.04722 0.02809 1.681 0.10162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06113 on 35 degrees of freedom
## Multiple R-squared: 0.4353, Adjusted R-squared: 0.3869
## F-statistic: 8.994 on 3 and 35 DF, p-value: 0.0001491
## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 2.63947 1 706.3882 < 2.2e-16 ***
## Site 0.10082 3 8.9943 0.0001491 ***
## Residuals 0.13078 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cover_rain.mod.jac)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.06192816 -0.13767379 0.013817463 0.1417750
## PFCA 2-TP -0.07925589 -0.15500152 -0.003510265 0.0374144
## PFCA 3-TP 0.04721930 -0.02852633 0.122964926 0.3485548
## PFCA 2-PFCA 1 -0.01732773 -0.09105311 0.056397651 0.9203909
## PFCA 3-PFCA 1 0.10914746 0.03542208 0.182872842 0.0017366
## PFCA 3-PFCA 2 0.12647519 0.05274981 0.200200570 0.0002772
# Contrasts
standardize <- emmeans::emmeans(Cover_rain.mod.bray, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.1340 0.0272 35 -4.925 0.0001
## PFCA 2 - TP -0.0707 0.0272 35 -2.598 0.0371
## PFCA 3 - TP 0.0231 0.0272 35 0.849 0.7122
##
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(Cover_rain.mod.bray, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.514 0.0197 35 0.474 0.554
## PFCA 1 0.380 0.0187 35 0.342 0.418
## PFCA 2 0.443 0.0187 35 0.405 0.481
## PFCA 3 0.537 0.0187 35 0.499 0.575
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.380, 0.443, 0.537, 0.514)
lower <- c(0.342, 0.405, 0.499, 0.474)
upper <- c(0.418, 0.481, 0.575, 0.554)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_rain_dissimilarity <- data.frame(site, response, lower, upper)
## 95% confidence intervals
# Contrasts
standardize <- emmeans::emmeans(Cover_rain.mod.jac, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.0619 0.0281 35 -2.205 0.0892
## PFCA 2 - TP -0.0793 0.0281 35 -2.822 0.0216
## PFCA 3 - TP 0.0472 0.0281 35 1.681 0.2427
##
## P value adjustment: dunnettx method for 3 tests
emmeans::emmeans(Cover_rain.mod.jac, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.542 0.0204 35 0.500 0.583
## PFCA 1 0.480 0.0193 35 0.440 0.519
## PFCA 2 0.462 0.0193 35 0.423 0.502
## PFCA 3 0.589 0.0193 35 0.550 0.628
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.480, 0.462, 0.589, 0.542)
lower <- c(0.440, .423, 0.550, 0.500)
upper <- c(0.519, 0.502, 0.628, 0.583)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_rain_jac_dissimilarity <- data.frame(site, response, lower, upper)
cover_rain_dis_wide <- cover_rain_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0))
## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.
cover_rain_dis_long <- cover_rain_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:177) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
cover_rain_dis_long_unique <- cover_rain_dis_long[!(cover_rain_dis_long$Aboveground == 0 & cover_rain_dis_long$Seed_Rain == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_rain_dis_long_unique <- cover_rain_dis_long_unique %>%
mutate(Aboveground = if_else(Aboveground > 0, 1, 0)) %>%
mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))
# Make a new column that adds both the seed rain and aboveground presence/absence columns. Then remove rows that = 2 (shares both species)
cover_rain_dis_long_unique$shared <- cover_rain_dis_long_unique$Aboveground + cover_rain_dis_long_unique$Seed_Rain
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique[!(cover_rain_dis_long_unique$shared == 2),]
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the seed rain
Unique_species_aboveground_com_rain <- cover_rain_dis_long_unique2[, -5]
Unique_species_aboveground_com_rain <- Unique_species_aboveground_com_rain %>%
filter(Aboveground != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Aboveground = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count unique taxa per site and transect for the seed rain compared to the seed bank
Unique_species_rain_com_aboveground <- cover_rain_dis_long_unique2[, -4]
Unique_species_rain_com_aboveground <- Unique_species_rain_com_aboveground %>%
filter(Seed_Rain != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Rain = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count the number of shared taxa between communities
cover_rain_dis_long_shared <- cover_rain_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
## Merge them together
Species_counts_cover_rain <- full_join(Unique_species_aboveground_com_rain, Unique_species_rain_com_aboveground) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
Species_counts_cover_rain <- full_join(Species_counts_cover_rain, cover_rain_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
Species_counts_cover_rain$Tot_Richness <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared
# Count the number of species in the seed rain for each transect
Species_counts_cover_rain$Tot_Rain <- Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared
# Count the number of species in the seed bank for each transect
Species_counts_cover_rain$Tot_Aboveground <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$Shared
## Immigrants
### What proportion of species dispersing in the seed rain are not found in the local cover?
Species_counts_cover_rain$Prop_Immigrant <- Species_counts_cover_rain$New_Rain / Species_counts_cover_rain$Tot_Rain
favstats(Species_counts_cover_rain$Prop_Immigrant~Species_counts_cover_rain$Site)
## Species_counts_cover_rain$Site min Q1 median Q3
## 1 PFCA 1 0.2500000 0.2660970 0.2891738 0.3793860
## 2 PFCA 2 0.2000000 0.2572674 0.2905844 0.3098639
## 3 PFCA 3 0.1818182 0.2826705 0.3675214 0.3942669
## 4 TP 0.2173913 0.2666667 0.3225806 0.3750000
## max mean sd n missing
## 1 0.5000000 0.3247644 0.08347182 10 0
## 2 0.3333333 0.2826759 0.04123406 10 0
## 3 0.4736842 0.3372457 0.09468692 10 0
## 4 0.4444444 0.3215489 0.07336399 9 0
# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?
# Species_counts_cover_rain.binom.mod <- glm(cbind(New_Rain, Tot_Rain) ~ Site, family = binomial, data =Species_counts_cover_rain)
# summary(Species_counts_cover_rain.binom.mod )
# Anova(Species_counts_cover_rain.binom.mod , type = "III")
Species_counts_cover_rain$Site <- factor(Species_counts_cover_rain$Site, levels=c("TP", "PFCA 2", "PFCA 3", "PFCA 1"))
Species_counts_cover_rain.lin.mod <- lm( Prop_Immigrant~ Site, data =Species_counts_cover_rain)
summary(Species_counts_cover_rain.lin.mod)
##
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_rain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.155428 -0.054688 0.003038 0.052547 0.175236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.321549 0.025308 12.705 1.14e-14 ***
## SitePFCA 2 -0.038873 0.034885 -1.114 0.273
## SitePFCA 3 0.015697 0.034885 0.450 0.656
## SitePFCA 1 0.003215 0.034885 0.092 0.927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07592 on 35 degrees of freedom
## Multiple R-squared: 0.07626, Adjusted R-squared: -0.002916
## F-statistic: 0.9632 on 3 and 35 DF, p-value: 0.421
Anova(Species_counts_cover_rain.lin.mod, type = "III")
## Anova Table (Type III tests)
##
## Response: Prop_Immigrant
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.93054 1 161.4255 1.141e-14 ***
## Site 0.01666 3 0.9632 0.421
## Residuals 0.20176 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Contrasts
standardize <- emmeans::emmeans(Species_counts_cover_rain.lin.mod, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 2 - TP -0.03887 0.0349 35 -1.114 0.5463
## PFCA 3 - TP 0.01570 0.0349 35 0.450 0.9162
## PFCA 1 - TP 0.00322 0.0349 35 0.092 0.9968
##
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(Species_counts_cover_rain.lin.mod, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.322 0.0253 35 0.270 0.373
## PFCA 2 0.283 0.0240 35 0.234 0.331
## PFCA 3 0.337 0.0240 35 0.289 0.386
## PFCA 1 0.325 0.0240 35 0.276 0.374
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.325, 0.283, 0.337, 0.322)
lower <- c(0.276, 0.234, 0.289, 0.270)
upper <- c( 0.374, 0.331, 0.386, 0.373)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Immigration_cover_rain <- data.frame(site, response, lower, upper)
#Re-order the levels
Species_counts_cover_rain$Site <- factor(Species_counts_cover_rain$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))
prop_new_species_cover_rain.plot <- ggplot()+
geom_jitter(data=Species_counts_cover_rain, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
geom_point(data = Immigration_cover_rain, aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
geom_errorbar(data = Immigration_cover_rain, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1)+
theme_classic()+
labs(x = "", y = "Prop. new species in seed rain")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(0, 1.075))+
geom_vline(xintercept = 3.5, linetype = "longdash")
prop_new_species_cover_rain.plot <- prop_new_species_cover_rain.plot + draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
x = .5, y = .5, width = .6)+
draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/SPHOBT_Sil.png",
x = 1.23, y = .5, width = .6) +
annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4)+
annotate("text", x = 4 , y = .97, label = "n.s.", size = 4.5)
prop_new_species_cover_rain.plot
names(cover_rain_dis_long) <- c("Site", "Transect", "SPP6", "Aboveground_Abundance", "Seed_Rain_Abundance")
cover_rain_dis_long_unique_abundance <- full_join(cover_rain_dis_long_unique, cover_rain_dis_long)
## Joining with `by = join_by(Site, Transect, SPP6)`
cover_rain_dis_long_shared_abundance <- cover_rain_dis_long_unique_abundance %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared_Seeds = sum(Seed_Rain_Abundance))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
cover_rain_dis_long_immigrant_abundance <- cover_rain_dis_long_unique_abundance %>%
filter(shared == 1) %>%
group_by(Site, Transect) %>%
summarize(Immigrant_Seeds = sum(Seed_Rain_Abundance))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
Prop_Immigrant_Seeds <- full_join(cover_rain_dis_long_shared_abundance, cover_rain_dis_long_immigrant_abundance)
## Joining with `by = join_by(Site, Transect)`
Prop_Immigrant_Seeds$Total_Seeds <- Prop_Immigrant_Seeds$Shared_Seeds + Prop_Immigrant_Seeds$Immigrant_Seeds
Prop_Immigrant_Seeds$Prop <- Prop_Immigrant_Seeds$Immigrant_Seeds / Prop_Immigrant_Seeds$Total_Seeds
Prop_Immigrant_Seeds$Site <- factor(Prop_Immigrant_Seeds$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))
# Prop_immigrant_cover_rain.binom.mod <- glm(cbind(Immigrant_Seeds, Total_Seeds) ~ Site, family = binomial, data =Prop_Immigrant_Seeds)
# summary(Prop_immigrant_cover_rain.binom.mod)
# Anova(Prop_immigrant_cover_rain.binom.mod, type = "III")
Prop_immigrant_cover_rain.lm.mod <- lm(Prop ~ Site, data = Prop_Immigrant_Seeds)
summary(Prop_immigrant_cover_rain.lm.mod)
##
## Call:
## lm(formula = Prop ~ Site, data = Prop_Immigrant_Seeds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.066032 -0.020370 -0.010473 0.005199 0.256420
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07927 0.01874 4.230 0.00016 ***
## SitePFCA 1 -0.06559 0.02583 -2.539 0.01571 *
## SitePFCA 2 -0.04865 0.02583 -1.884 0.06795 .
## SitePFCA 3 -0.04254 0.02583 -1.647 0.10855
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05622 on 35 degrees of freedom
## Multiple R-squared: 0.1644, Adjusted R-squared: 0.09274
## F-statistic: 2.295 on 3 and 35 DF, p-value: 0.09484
Anova(Prop_immigrant_cover_rain.lm.mod, type = "3")
## Anova Table (Type III tests)
##
## Response: Prop
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.056559 1 17.8957 0.0001595 ***
## Site 0.021758 3 2.2948 0.0948394 .
## Residuals 0.110616 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Contrasts
standardize <- emmeans::emmeans(Prop_immigrant_cover_rain.lm.mod, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.0656 0.0258 35 -2.539 0.0426
## PFCA 2 - TP -0.0487 0.0258 35 -1.884 0.1691
## PFCA 3 - TP -0.0425 0.0258 35 -1.647 0.2572
##
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(Prop_immigrant_cover_rain.lm.mod, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.0793 0.0187 35 0.041231 0.1173
## PFCA 1 0.0137 0.0178 35 -0.022403 0.0498
## PFCA 2 0.0306 0.0178 35 -0.005471 0.0667
## PFCA 3 0.0367 0.0178 35 0.000645 0.0728
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.0137, 0.0306, 0.0367, 0.0793)
lower <- c(-0.022403, -0.005471, 0.000645, 0.041231)
upper <- c(0.0498, 0.0667, 0.0728, 0.1173)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Prop_Immigration_cover_rain <- data.frame(site, response, lower, upper)
#Re-order the levels
Prop_Immigrant_Seeds$Site <- factor(Prop_Immigrant_Seeds$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))
prop_immigrant_seeds_cover_rain.plot <- ggplot()+
geom_jitter(data=Prop_Immigrant_Seeds, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
theme_classic()+
geom_point(data = Prop_Immigration_cover_rain, aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
geom_errorbar(data = Prop_Immigration_cover_rain, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1)+
labs(x = "", y = "Prop. new seeds in seed rain", title = "")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(-0.025, 1.075))+
geom_vline(xintercept = 3.5, linetype = "longdash")
# The two "outlier" points are because we caught a lot of Solidago altissima and Barbarea vulgaris at these points but they weren't in the immediate cover.
prop_immigrant_seeds_cover_rain.plot <- prop_immigrant_seeds_cover_rain.plot +
draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
x = .5, y = .5, width = .6)+
draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/SPHOBT_Sil.png",
x = 1.23, y = .5, width = .6)+
annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4) +
annotate("text", x = 1 , y = .1, label = "*", size = 9.5)
prop_immigrant_seeds_cover_rain.plot
cover_bank_sum <- cover_bank %>%
group_by(SPP6) %>%
summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seedlings))
Germinated but didn’t observed in the cover
Almost everything we were able to germinate that was not in the local cover were weedy annual species. Some things like the Cardamine might have been in the cover but have such an early phenology that we missed them. The two tree species might be greenhouse contaminates since the roof was partially open and these trees were right next door.
Observed in the cover but didn’t germinate
Too many to list (112 taxa)
Clearly, only a small fraction of local species survive the transition from cover -> seed rain -> to germinable seed bank.
Few perennial species were present in the germinable seed bank.
## # A tibble: 9,984 × 5
## # Groups: Site, Transect, SPP6 [7,605]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ABUTHE Seed_Bank 0
## 2 PFCA 1 1 ACAVIR Aboveground 11.5
## 3 PFCA 1 1 ACAVIR Seed_Bank 0
## 4 PFCA 1 1 ACERUB Aboveground 0
## 5 PFCA 1 1 ACHMIL Aboveground 0
## 6 PFCA 1 1 ACHMIL Seed_Bank 0
## 7 PFCA 1 1 AGAFAS Aboveground 0
## 8 PFCA 1 1 AGATEN Aboveground 0
## 9 PFCA 1 1 AGRGIG Aboveground 0
## 10 PFCA 1 1 AGRHYE Aboveground 0
## # ℹ 9,974 more rows
# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two
##### run dissimilarity loops
cover_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(cover_bank_dis$Site))
for(s in 1:length(sites)){
temp <-subset(cover_bank_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
#for(b in 1:length(blocks)){
#temp_b<-subset(temp, block==blocks[b])
#plots<-as.vector(unique(temp_b$plot))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"hellinger"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
cover_bank_dis_results <-rbind( cover_bank_dis_results, new.row)
}
print(s/length(sites))
}
## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1
row.names(cover_bank_dis_results )<-NULL
cover_bank_dis_results <- data.frame(cover_bank_dis_results )
names(cover_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")
##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = cover_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.187831 -0.029607 -0.002413 0.043986 0.171576
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.765696 0.023016 33.268 < 2e-16 ***
## SitePFCA 1 -0.192711 0.031725 -6.074 6.16e-07 ***
## SitePFCA 2 -0.007888 0.031725 -0.249 0.805
## SitePFCA 3 0.044227 0.031725 1.394 0.172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06905 on 35 degrees of freedom
## Multiple R-squared: 0.6634, Adjusted R-squared: 0.6345
## F-statistic: 22.99 on 3 and 35 DF, p-value: 2.111e-08
## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 5.2766 1 1106.754 < 2.2e-16 ***
## Site 0.3289 3 22.993 2.111e-08 ***
## Residuals 0.1669 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cover_bank.mod.bray)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.192711214 -0.27827161 -0.10715081 0.0000036
## PFCA 2-TP -0.007887543 -0.09344794 0.07767286 0.9945027
## PFCA 3-TP 0.044227483 -0.04133292 0.12978788 0.5113778
## PFCA 2-PFCA 1 0.184823671 0.10154529 0.26810205 0.0000047
## PFCA 3-PFCA 1 0.236938697 0.15366032 0.32021707 0.0000000
## PFCA 3-PFCA 2 0.052115027 -0.03116335 0.13539340 0.3452273
##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = cover_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.108134 -0.051282 -0.007532 0.055907 0.133446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.83530 0.02209 37.818 <2e-16 ***
## SitePFCA 1 -0.05277 0.03045 -1.733 0.0918 .
## SitePFCA 2 -0.00717 0.03045 -0.236 0.8152
## SitePFCA 3 0.01687 0.03045 0.554 0.5829
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06626 on 35 degrees of freedom
## Multiple R-squared: 0.1469, Adjusted R-squared: 0.07374
## F-statistic: 2.008 on 3 and 35 DF, p-value: 0.1307
## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.2796 1 1430.1902 <2e-16 ***
## Site 0.0265 3 2.0083 0.1307
## Residuals 0.1537 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = x)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.052771996 -0.13488081 0.02933682 0.3222879
## PFCA 2-TP -0.007169957 -0.08927877 0.07493886 0.9953158
## PFCA 3-TP 0.016874135 -0.06523468 0.09898295 0.9447732
## PFCA 2-PFCA 1 0.045602039 -0.03431681 0.12552089 0.4259519
## PFCA 3-PFCA 1 0.069646131 -0.01027272 0.14956498 0.1060862
## PFCA 3-PFCA 2 0.024044092 -0.05587476 0.10396294 0.8486613
# Contrasts
standardize <- emmeans::emmeans(Cover_bank.mod.bray, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.19271 0.0317 35 -6.074 <.0001
## PFCA 2 - TP -0.00789 0.0317 35 -0.249 0.9752
## PFCA 3 - TP 0.04423 0.0317 35 1.394 0.3808
##
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(Cover_bank.mod.bray, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.766 0.0230 35 0.719 0.812
## PFCA 1 0.573 0.0218 35 0.529 0.617
## PFCA 2 0.758 0.0218 35 0.713 0.802
## PFCA 3 0.810 0.0218 35 0.766 0.854
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.573, 0.758, 0.810, 0.766)
lower <- c(0.529, 0.714, 0.766, 0.719)
upper <- c(0.617, 0.802, 0.854, 0.812)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_bank_dissimilarity <- data.frame(site, response, lower, upper)
## 95% confidence intervals
standardize <- emmeans::emmeans(Cover_bank.mod.jac, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.05277 0.0304 35 -1.733 0.2218
## PFCA 2 - TP -0.00717 0.0304 35 -0.236 0.9778
## PFCA 3 - TP 0.01687 0.0304 35 0.554 0.8725
##
## P value adjustment: dunnettx method for 3 tests
emmeans::emmeans(Cover_bank.mod.jac, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.835 0.0221 35 0.790 0.880
## PFCA 1 0.783 0.0210 35 0.740 0.825
## PFCA 2 0.828 0.0210 35 0.786 0.871
## PFCA 3 0.852 0.0210 35 0.810 0.895
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.783, 0.828, 0.852, 0.835)
lower <- c(0.740, .786, 0.810, 0.790)
upper <- c( 0.825, 0.871, 0.895, 0.880)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)
cover_bank_dis_wide <- cover_bank_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0))
## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.
cover_bank_dis_long <- cover_bank_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:198) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
cover_bank_dis_long_unique <- cover_bank_dis_long[!(cover_bank_dis_long$Seed_Bank == 0 & cover_bank_dis_long$Aboveground == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_bank_dis_long_unique <- cover_bank_dis_long_unique %>%
mutate(Seed_Bank = if_else(Seed_Bank > 0, 1, 0)) %>%
mutate(Aboveground = if_else(Aboveground > 0, 1, 0))
# Make a new column that adds both the cover and seed bank presence/absence columns. Then remove rows that = 2 (shares both species)
cover_bank_dis_long_unique$shared <- cover_bank_dis_long_unique$Seed_Bank + cover_bank_dis_long_unique$Aboveground
cover_bank_dis_long_unique2 <- cover_bank_dis_long_unique[!(cover_bank_dis_long_unique$shared == 2),]
cover_bank_dis_long_unique2 <- cover_bank_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the cover
Unique_species_bank_com_cover <- cover_bank_dis_long_unique2[, -4]
Unique_species_bank_com_cover <- Unique_species_bank_com_cover %>%
filter(Seed_Bank != 0 ) %>% group_by(Site, Transect) %>%
summarize(New_Bank = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count unique taxa per site and transect for the cover compared to the seed bank
Unique_species_cover_com_bank <- cover_bank_dis_long_unique2[, -5]
Unique_species_cover_com_bank <- Unique_species_cover_com_bank %>%
filter(Aboveground != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Cover = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
# Count the number of shared taxa between communities
cover_bank_dis_long_shared <- cover_bank_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
## Merge them together
Species_counts_cover_bank <- full_join(Unique_species_bank_com_cover, Unique_species_cover_com_bank) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
Species_counts_cover_bank <- full_join(Species_counts_cover_bank, cover_bank_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0))
## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
Species_counts_cover_bank$Tot_Richness <- Species_counts_cover_bank$New_Bank + Species_counts_cover_bank$New_Cover + Species_counts_cover_bank$Shared
# Count the number of species in the seed rain for each transect
Species_counts_cover_bank$Tot_Cover <- Species_counts_cover_bank$New_Cover + Species_counts_cover_bank$Shared
# Count the number of species in the seed bank for each transect
Species_counts_cover_bank$Tot_Bank <- Species_counts_cover_bank$New_Bank + Species_counts_cover_bank$Shared
## Immigrants
Species_counts_cover_bank$Prop_Immigrant <- Species_counts_cover_bank$New_Bank / Species_counts_cover_bank$Tot_Bank
# No difference in proportion of immigrant species
#Re-order the levels
Species_counts_cover_bank$Site <- factor(Species_counts_cover_bank$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))
Species_counts_cover_bank.lin.mod <- lm( Prop_Immigrant~ Site, data = Species_counts_cover_bank)
summary(Species_counts_cover_bank.lin.mod)
##
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_bank)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.266163 -0.091067 -0.006404 0.084063 0.283837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.46616 0.04369 10.669 1.51e-12 ***
## SitePFCA 1 -0.12643 0.06022 -2.099 0.0431 *
## SitePFCA 2 -0.05650 0.06022 -0.938 0.3546
## SitePFCA 3 0.10631 0.06022 1.765 0.0863 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1311 on 35 degrees of freedom
## Multiple R-squared: 0.3251, Adjusted R-squared: 0.2673
## F-statistic: 5.621 on 3 and 35 DF, p-value: 0.002972
Anova(Species_counts_cover_bank.lin.mod, type = "III")
## Anova Table (Type III tests)
##
## Response: Prop_Immigrant
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.95578 1 113.8357 1.515e-12 ***
## Site 0.28971 3 5.6209 0.002972 **
## Residuals 0.60132 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Contrasts
standardize <- emmeans::emmeans(Species_counts_cover_bank.lin.mod, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.1264 0.0602 35 -2.099 0.1110
## PFCA 2 - TP -0.0565 0.0602 35 -0.938 0.6572
## PFCA 3 - TP 0.1063 0.0602 35 1.765 0.2098
##
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(Species_counts_cover_bank.lin.mod, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.466 0.0437 35 0.377 0.555
## PFCA 1 0.340 0.0414 35 0.256 0.424
## PFCA 2 0.410 0.0414 35 0.326 0.494
## PFCA 3 0.572 0.0414 35 0.488 0.657
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.340, 0.410, 0.572, 0.466)
lower <- c(0.256, 0.326, 0.488 , 0.377)
upper <- c( 0.424, 0.494, 0.657, 0.555)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Immigration_cover_bank <- data.frame(site, response, lower, upper)
### the oldest restoration has the greatest proportion of immigrants -> likely caused by a persistent seed bank formed from when this site was farmed prior to restoration. Weedy non-native species likely were displaced from the aboveground vegetation but remained present in the seed bank. This trend is not seen in Tucker Prairie which doesn't have the same past.
Species_counts_cover_bank$Site <- factor(Species_counts_cover_bank$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))
prop_new_species_cover_bank.plot <- ggplot() +
geom_jitter(data = Species_counts_cover_bank, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = 0.2, size = 2.5) +
geom_point(data = Immigration_cover_bank , aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
geom_errorbar(data = Immigration_cover_bank , aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1) +
theme_classic()+
labs(x = "", y = "Prop. new species in seed bank")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(0, 1.075))+
geom_vline(xintercept = 3.5, linetype = "longdash")
prop_new_species_cover_bank.plot <- prop_new_species_cover_bank.plot + draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
x = .5, y = .5, width = .6)+
draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/AMBART_Sil.png",
x = 1.15, y = .5, width = .75) +
annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4)+
annotate("text", x = 4 , y = .97, label = "n.s.", size = 4.5)
prop_new_species_cover_bank.plot
names(cover_bank_dis_long) <- c("Site", "Transect", "SPP6", "Aboveground_Abundance", "Seed_Bank_Abundance")
cover_bank_dis_long_unique_abundance <- full_join(cover_bank_dis_long_unique, cover_bank_dis_long)
## Joining with `by = join_by(Site, Transect, SPP6)`
cover_bank_dis_long_shared_abundance <- cover_bank_dis_long_unique_abundance %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared_Seedlings = sum(Seed_Bank_Abundance))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
cover_bank_dis_long_immigrant_abundance <- cover_bank_dis_long_unique_abundance %>%
filter(shared == 1) %>%
group_by(Site, Transect) %>%
summarize(Immigrant_Seedlings = sum(Seed_Bank_Abundance))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
Prop_Immigrant_Seedlings <- full_join(cover_bank_dis_long_shared_abundance, cover_bank_dis_long_immigrant_abundance)
## Joining with `by = join_by(Site, Transect)`
Prop_Immigrant_Seedlings$Total_Seedlings <- Prop_Immigrant_Seedlings$Shared_Seedlings + Prop_Immigrant_Seedlings$Immigrant_Seedlings
Prop_Immigrant_Seedlings$Prop <- Prop_Immigrant_Seedlings$Immigrant_Seedlings / Prop_Immigrant_Seedlings$Total_Seedlings
# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?
## binomial regression produces results that don't work with the actual observed data. Linear model works just as well with more appropriate mean and CI results.
#Re-order the levels
Prop_Immigrant_Seedlings$Site <- factor(Prop_Immigrant_Seedlings$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))
Prop_Immigrant_Seedlings.lin.mod <- lm(Prop ~ Site, data = Prop_Immigrant_Seedlings)
summary(Prop_Immigrant_Seedlings.lin.mod)
##
## Call:
## lm(formula = Prop ~ Site, data = Prop_Immigrant_Seedlings)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33648 -0.14610 -0.00805 0.08538 0.34960
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.45040 0.05765 7.812 3.53e-09 ***
## SitePFCA 1 -0.39907 0.07947 -5.022 1.50e-05 ***
## SitePFCA 2 -0.19101 0.07947 -2.404 0.0217 *
## SitePFCA 3 0.11288 0.07947 1.420 0.1643
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.173 on 35 degrees of freedom
## Multiple R-squared: 0.5891, Adjusted R-squared: 0.5539
## F-statistic: 16.72 on 3 and 35 DF, p-value: 6.561e-07
Anova(Prop_Immigrant_Seedlings.lin.mod, type = "III")
## Anova Table (Type III tests)
##
## Response: Prop
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.8258 1 61.032 3.529e-09 ***
## Site 1.5009 3 16.725 6.561e-07 ***
## Residuals 1.0470 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Contrasts
standardize <- emmeans::emmeans(Prop_Immigrant_Seedlings.lin.mod , "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.399 0.0795 35 -5.022 <.0001
## PFCA 2 - TP -0.191 0.0795 35 -2.404 0.0579
## PFCA 3 - TP 0.113 0.0795 35 1.420 0.3666
##
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(Prop_Immigrant_Seedlings.lin.mod, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.4504 0.0577 35 0.3334 0.567
## PFCA 1 0.0513 0.0547 35 -0.0597 0.162
## PFCA 2 0.2594 0.0547 35 0.1484 0.370
## PFCA 3 0.5633 0.0547 35 0.4522 0.674
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.0513, 0.2594, 0.5633, 0.4504)
lower <- c(-0.0597, 0.1484, 0.4522 , 0.3334)
upper <- c( 0.162, 0.370, 0.674, 0.567)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Prop_cover_bank <- data.frame(site, response, lower, upper)
#Re-order the levels
Prop_Immigrant_Seedlings$Site <- factor(Prop_Immigrant_Seedlings$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))
prop_immigrant_seedlings_cover_bank.plot <- ggplot()+
geom_jitter(data=Prop_Immigrant_Seedlings, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
geom_point(data = Prop_cover_bank, aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
geom_errorbar(data = Prop_cover_bank, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1) +
theme_classic()+
labs(x = "", y = "Prop. new seeds in seed bank", title = "")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(-0.06, 1.075))+
geom_vline(xintercept = 3.5, linetype = "longdash")
prop_immigrant_seedlings_cover_bank.plot <- prop_immigrant_seedlings_cover_bank.plot + draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
x = .5, y = .5, width = .6)+
draw_image(
"~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/AMBART_Sil.png",
x = 1.15, y = .5, width = .75) +
annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4)+
annotate("text", x = 1 , y = .4, label = "*", size = 9.5)
rain_bank_sum <- rain_bank %>%
group_by(SPP6) %>%
summarize(tot.seeds = sum(Number_Seeds), tot.seedlings = sum(Number_Seedlings))
Germinated but didn’t catch in the seed rain
Mostly ruderal annual species.
Caught in the seed rain but didn’t germinate
Too many to list (61 taxa)
## # A tibble: 9,984 × 5
## # Groups: Site, Transect, SPP6 [7,605]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ABUTHE Seed_Bank 0
## 2 PFCA 1 1 ACAVIR Aboveground 11.5
## 3 PFCA 1 1 ACAVIR Seed_Bank 0
## 4 PFCA 1 1 ACERUB Aboveground 0
## 5 PFCA 1 1 ACHMIL Aboveground 0
## 6 PFCA 1 1 ACHMIL Seed_Bank 0
## 7 PFCA 1 1 AGAFAS Aboveground 0
## 8 PFCA 1 1 AGATEN Aboveground 0
## 9 PFCA 1 1 AGRGIG Aboveground 0
## 10 PFCA 1 1 AGRHYE Aboveground 0
## # ℹ 9,974 more rows
## # A tibble: 8,200 × 5
## # Groups: Site, Transect, SPP6 [5,600]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ABUTHE Seed_Bank 0
## 2 PFCA 1 1 ACAVIR Seed_Bank 0
## 3 PFCA 1 1 ACAVIR Seed_Rain 22
## 4 PFCA 1 1 ACHMIL Seed_Bank 0
## 5 PFCA 1 1 ACHMIL Seed_Rain 0
## 6 PFCA 1 1 AGASPP Seed_Rain 0
## 7 PFCA 1 1 AGRHYE Seed_Bank 0
## 8 PFCA 1 1 AGRHYE Seed_Rain 0
## 9 PFCA 1 1 ALOCAR Seed_Bank 0
## 10 PFCA 1 1 ALOCAR Seed_Rain 0
## # ℹ 8,190 more rows
# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two
##### run dissimilarity loops
rain_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(rain_bank_dis$Site))
for(s in 1:length(sites)){
temp <-subset(rain_bank_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"hellinger"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
rain_bank_dis_results <-rbind(rain_bank_dis_results, new.row)
}
print(s/length(sites))
}
## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1
row.names(rain_bank_dis_results)<-NULL
rain_bank_dis_results <- data.frame(rain_bank_dis_results)
names(rain_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")
##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = rain_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.170306 -0.047042 -0.002914 0.046581 0.125648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.61439 0.02451 25.063 <2e-16 ***
## SitePFCA 1 -0.11376 0.03467 -3.281 0.0023 **
## SitePFCA 2 0.05321 0.03467 1.535 0.1336
## SitePFCA 3 0.06543 0.03467 1.887 0.0672 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07752 on 36 degrees of freedom
## Multiple R-squared: 0.481, Adjusted R-squared: 0.4377
## F-statistic: 11.12 on 3 and 36 DF, p-value: 2.602e-05
## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 3.7748 1 628.18 < 2.2e-16 ***
## Site 0.2005 3 11.12 2.602e-05 ***
## Residuals 0.2163 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rain_bank.mod.bray)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.11375660 -0.20712329 -0.02038992 0.0117875
## PFCA 2-TP 0.05320739 -0.04015930 0.14657407 0.4279537
## PFCA 3-TP 0.06543023 -0.02793646 0.15879691 0.2513341
## PFCA 2-PFCA 1 0.16696399 0.07359731 0.26033068 0.0001492
## PFCA 3-PFCA 1 0.17918683 0.08582015 0.27255352 0.0000514
## PFCA 3-PFCA 2 0.01222284 -0.08114385 0.10558953 0.9847184
## 95% confidence intervals
# Contrasts
standardize <- emmeans::emmeans(Rain_bank.mod.bray, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.1138 0.0347 36 -3.281 0.0065
## PFCA 2 - TP 0.0532 0.0347 36 1.535 0.3079
## PFCA 3 - TP 0.0654 0.0347 36 1.887 0.1674
##
## P value adjustment: dunnettx method for 3 tests
emmeans::emmeans(Rain_bank.mod.bray, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.614 0.0245 36 0.565 0.664
## PFCA 1 0.501 0.0245 36 0.451 0.550
## PFCA 2 0.668 0.0245 36 0.618 0.717
## PFCA 3 0.680 0.0245 36 0.630 0.730
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.501, 0.668, 0.680, 0.614)
lower <- c(0.451, 0.618, 0.630, 0.565)
upper <- c(0.550, 0.717, 0.730, 0.664)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Rain_bank_dissimilarity <- data.frame(site, response, lower, upper)
#Re-order the levels
rain_bank_dis_results$Site <- factor(rain_bank_dis_results$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))
Rain_bank.mod.jac <- lm(dissimilarity_jaccard ~ Site, data= rain_bank_dis_results)
summary(Rain_bank.mod.jac)
##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = rain_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16845 -0.03871 0.00568 0.02850 0.13258
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.76505 0.02153 35.540 <2e-16 ***
## SitePFCA 1 -0.04397 0.03044 -1.444 0.157
## SitePFCA 2 -0.01655 0.03044 -0.544 0.590
## SitePFCA 3 0.04111 0.03044 1.350 0.185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06807 on 36 degrees of freedom
## Multiple R-squared: 0.1856, Adjusted R-squared: 0.1178
## F-statistic: 2.735 on 3 and 36 DF, p-value: 0.05776
Anova(Rain_bank.mod.jac, type = "III")
## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 5.8529 1 1263.0564 < 2e-16 ***
## Site 0.0380 3 2.7354 0.05776 .
## Residuals 0.1668 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 95% confidence intervals
# Contrasts
standardize <- emmeans::emmeans(Rain_bank.mod.jac, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.0440 0.0304 36 -1.444 0.3536
## PFCA 2 - TP -0.0166 0.0304 36 -0.544 0.8772
## PFCA 3 - TP 0.0411 0.0304 36 1.350 0.4047
##
## P value adjustment: dunnettx method for 3 tests
emmeans::emmeans(Rain_bank.mod.jac, "Site", type = "response")
## Site emmean SE df lower.CL upper.CL
## TP 0.765 0.0215 36 0.721 0.809
## PFCA 1 0.721 0.0215 36 0.677 0.765
## PFCA 2 0.748 0.0215 36 0.705 0.792
## PFCA 3 0.806 0.0215 36 0.762 0.850
##
## Confidence level used: 0.95
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.721, 0.748, 0.806, 0.765)
lower <- c(0.677, 0.705, 0.762, 0.721)
upper <- c(0.765, 0.792, 0.850, 0.809)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Rain_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)
## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
## There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
##
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_rain_bank)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.288874 -0.079552 -0.009654 0.067853 0.220249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28887 0.03892 7.423 9.23e-09 ***
## SitePFCA 1 -0.07665 0.05504 -1.393 0.172269
## SitePFCA 2 -0.06449 0.05504 -1.172 0.248973
## SitePFCA 3 0.22421 0.05504 4.074 0.000243 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1231 on 36 degrees of freedom
## Multiple R-squared: 0.5179, Adjusted R-squared: 0.4778
## F-statistic: 12.89 on 3 and 36 DF, p-value: 7.119e-06
## Anova Table (Type III tests)
##
## Response: Prop_Immigrant
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.83448 1 55.098 9.227e-09 ***
## Site 0.58579 3 12.893 7.119e-06 ***
## Residuals 0.54524 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.0766 0.055 36 -1.393 0.3812
## PFCA 2 - TP -0.0645 0.055 36 -1.172 0.5104
## PFCA 3 - TP 0.2242 0.055 36 4.074 0.0007
##
## P value adjustment: dunnettx method for 3 tests
## Site emmean SE df lower.CL upper.CL
## TP 0.289 0.0389 36 0.210 0.368
## PFCA 1 0.212 0.0389 36 0.133 0.291
## PFCA 2 0.224 0.0389 36 0.145 0.303
## PFCA 3 0.513 0.0389 36 0.434 0.592
##
## Confidence level used: 0.95
##
## Call:
## lm(formula = Prop ~ Site, data = Prop_Immigrant_Seedlings_Rain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32186 -0.07339 -0.02271 0.03860 0.47470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28530 0.05768 4.946 1.77e-05 ***
## SitePFCA 1 -0.25212 0.08157 -3.091 0.00384 **
## SitePFCA 2 -0.18605 0.08157 -2.281 0.02857 *
## SitePFCA 3 0.17718 0.08157 2.172 0.03651 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1824 on 36 degrees of freedom
## Multiple R-squared: 0.4845, Adjusted R-squared: 0.4415
## F-statistic: 11.28 on 3 and 36 DF, p-value: 2.312e-05
## Anova Table (Type III tests)
##
## Response: Prop
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.81395 1 24.466 1.770e-05 ***
## Site 1.12545 3 11.277 2.312e-05 ***
## Residuals 1.19766 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast estimate SE df t.ratio p.value
## PFCA 1 - TP -0.252 0.0816 36 -3.091 0.0108
## PFCA 2 - TP -0.186 0.0816 36 -2.281 0.0754
## PFCA 3 - TP 0.177 0.0816 36 2.172 0.0951
##
## P value adjustment: dunnettx method for 3 tests
## Site emmean SE df lower.CL upper.CL
## TP 0.2853 0.0577 36 0.1683 0.402
## PFCA 1 0.0332 0.0577 36 -0.0838 0.150
## PFCA 2 0.0992 0.0577 36 -0.0177 0.216
## PFCA 3 0.4625 0.0577 36 0.3455 0.579
##
## Confidence level used: 0.95
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMATUB"
## [9] "AMBART" "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
## [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
## [25] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "COMUMB"
## [33] "CONCAN" "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPSPP" "DALPUR"
## [41] "DAUCAR" "DESILL" "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
## [49] "ELESPP" "ELYCAN" "ELYVIR" "ERASPE" "EREHIE" "ERISPP" "ERIVIL" "ERYYUC"
## [57] "EUPCOR" "EUPSPP" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT"
## [65] "GENSPP" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HORPUS" "HYPHIR"
## [73] "HYPPUN" "IPOHED" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPDEN"
## [81] "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LEUVUL" "LIAPYC" "LIASPP"
## [89] "LIASQU" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP"
## [97] "MOLVER" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENSPP" "OXADIL" "PANCAP"
## [105] "PANVIR" "PARINT" "PENDIG" "PERLON" "PLAMAJ" "PLAOCC" "PLAVIR" "POAPRA"
## [113] "POLSPP" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP"
## [153] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [161] "ULMPUM" "ULMSPP" "VERHAS" "VERHEL" "VEROSP" "VERSPP" "VIOSAG" "VIOSPP"
## [169] "VULOCT" "ZIZAUR" "PLAPUS" "POACHA" "ABUTHE" "POPDEL" "EUPHUM" "PANDIC"
## [177] "POROLE" "VERTHA" "DIPLAC" "DIGSAN" "TOXRAD"
### Number of seedlings at each site
Site_seedlings <- everything_merged %>%
filter(Type == "Seed_Bank") %>%
group_by(Site) %>%
summarize(tot.seedlings = sum(Abundance))
### Number of seedlings at each transect and site
Just_seedlings <- everything_merged %>%
filter(Type == "Seed_Bank") %>%
group_by(Site, Transect) %>%
summarize(tot.seedlings = sum(Abundance))
## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.
## Mean number and standard deviation of seedlings
Just_seedlings_mean <- Just_seedlings %>%
group_by(Site) %>%
summarize(mean.seedlings = mean(tot.seedlings), sd.seedlings = sd(tot.seedlings))
#Re-order the levels
Just_seedlings$Site <- factor(Just_seedlings$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))
# Fit a negative binomial GLM to predict number of seedlings as a function of site
seedlings.mod <- glm.nb(tot.seedlings ~ Site, data = Just_seedlings)
summary(seedlings.mod)
##
## Call:
## glm.nb(formula = tot.seedlings ~ Site, data = Just_seedlings,
## init.theta = 2.995739929, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.85486 0.19860 19.410 < 2e-16 ***
## SitePFCA 1 1.89803 0.27045 7.018 2.25e-12 ***
## SitePFCA 2 0.91920 0.27142 3.387 0.000707 ***
## SitePFCA 3 0.07303 0.27348 0.267 0.789434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.9957) family taken to be 1)
##
## Null deviance: 111.465 on 38 degrees of freedom
## Residual deviance: 40.926 on 35 degrees of freedom
## AIC: 428.26
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.996
## Std. Err.: 0.672
##
## 2 x log-likelihood: -418.256
Anova(seedlings.mod, type = "3")
## Analysis of Deviance Table (Type III tests)
##
## Response: tot.seedlings
## LR Chisq Df Pr(>Chisq)
## Site 70.539 3 3.272e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Contrasts
standardize <- emmeans::emmeans(seedlings.mod, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")
## contrast estimate SE df z.ratio p.value
## PFCA 1 - TP 1.898 0.270 Inf 7.018 <.0001
## PFCA 2 - TP 0.919 0.271 Inf 3.387 0.0021
## PFCA 3 - TP 0.073 0.273 Inf 0.267 0.9712
##
## Results are given on the log (not the response) scale.
## P value adjustment: dunnettx method for 3 tests
## 95% confidence intervals
emmeans::emmeans(seedlings.mod, "Site", type = "response")
## Site response SE df asymp.LCL asymp.UCL
## TP 47.2 9.38 Inf 32.0 69.7
## PFCA 1 315.1 57.84 Inf 219.9 451.5
## PFCA 2 118.4 21.90 Inf 82.4 170.1
## PFCA 3 50.8 9.55 Inf 35.1 73.4
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c( 315.1, 118.4, 50.8, 47.2)
lower <- c(219.9, 82.4, 35.1, 32.0)
upper <- c(451.5, 170.1, 73.4, 69.7)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
seedlings.mod.est <- data.frame(site, response, lower, upper)
#Re-order the levels
Just_seedlings$Site <- factor(Just_seedlings$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))
### Plot number of seedlings vs. site
ggplot(Just_seedlings, aes(x= Site, y = tot.seedlings, color = Site))+
geom_jitter(width =0.2, size = 2, show.legend = FALSE)+
geom_point(data = Just_seedlings_mean, aes(x = Site, y = mean.seedlings),size = 3.5, color = "black")+
geom_errorbar(data = seedlings.mod.est, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1, color = "black")+
theme_classic()+
labs(x= "Site", y = "Number of Seedlings")+
theme(text=element_text(size=18), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
scale_color_manual(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
labels = c("Young", "Middle", "Old", "Remnant"),
values = c("#1f78b4", "#a6cee3", "#b2df8a","#33a02c"))+
scale_x_discrete(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
labels = c("Young", "Middle", "Old", "Remnant"))+
geom_vline(xintercept = 3.5, linetype = "longdash")+
annotate("text", x = 1 , y = 850, label = "*", size = 9.5)+
annotate("text", x = 2 , y = 300, label = "*", size = 9.5)+
ylim(0, 1000)
## # A tibble: 39 × 3
## # Groups: Site [4]
## Site Transect Species_richness
## <fct> <fct> <int>
## 1 TP 1 34
## 2 TP 10 37
## 3 TP 2 30
## 4 TP 3 31
## 5 TP 4 27
## 6 TP 5 30
## 7 TP 6 29
## 8 TP 8 38
## 9 TP 9 43
## 10 PFCA 1 1 28
## # ℹ 29 more rows
##
## Call:
## glm(formula = Species_richness ~ Site, family = "poisson", data = Uni_Above_Diversity)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.503219 0.057831 60.576 < 2e-16 ***
## SitePFCA 1 0.121122 0.077532 1.562 0.118
## SitePFCA 2 0.314493 0.074447 4.224 2.4e-05 ***
## SitePFCA 3 0.002338 0.079671 0.029 0.977
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 66.379 on 38 degrees of freedom
## Residual deviance: 40.952 on 35 degrees of freedom
## AIC: 261.23
##
## Number of Fisher Scoring iterations: 4
## Analysis of Deviance Table (Type III tests)
##
## Response: Species_richness
## LR Chisq Df Pr(>Chisq)
## Site 25.427 3 1.257e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast estimate SE df z.ratio p.value
## PFCA 1 - TP 0.12112 0.0775 Inf 1.562 0.2796
## PFCA 2 - TP 0.31449 0.0744 Inf 4.224 0.0001
## PFCA 3 - TP 0.00234 0.0797 Inf 0.029 0.9997
##
## Results are given on the log (not the response) scale.
## P value adjustment: dunnettx method for 3 tests
## Site rate SE df asymp.LCL asymp.UCL
## TP 33.2 1.92 Inf 29.7 37.2
## PFCA 1 37.5 1.94 Inf 33.9 41.5
## PFCA 2 45.5 2.13 Inf 41.5 49.9
## PFCA 3 33.3 1.82 Inf 29.9 37.1
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## # A tibble: 39 × 3
## # Groups: Site [4]
## Site Transect Species_richness
## <fct> <fct> <int>
## 1 TP 1 34
## 2 TP 10 37
## 3 TP 2 30
## 4 TP 3 31
## 5 TP 4 27
## 6 TP 5 30
## 7 TP 6 29
## 8 TP 8 38
## 9 TP 9 43
## 10 PFCA 1 1 28
## # ℹ 29 more rows
##
## Call:
## glm(formula = Species_richness_native ~ Site, family = "poisson",
## data = Uni_Above_Diversity_native)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.44113 0.05965 57.684 < 2e-16 ***
## SitePFCA 1 -0.09123 0.08407 -1.085 0.27786
## SitePFCA 2 0.20953 0.07846 2.670 0.00757 **
## SitePFCA 3 -0.06014 0.08343 -0.721 0.47102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 58.191 on 38 degrees of freedom
## Residual deviance: 39.489 on 35 degrees of freedom
## AIC: 253.49
##
## Number of Fisher Scoring iterations: 4
## Analysis of Deviance Table (Type III tests)
##
## Response: Species_richness_native
## LR Chisq Df Pr(>Chisq)
## Site 18.703 3 0.000315 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast estimate SE df z.ratio p.value
## PFCA 1 - TP -0.0912 0.0841 Inf -1.085 0.5567
## PFCA 2 - TP 0.2095 0.0785 Inf 2.670 0.0213
## PFCA 3 - TP -0.0601 0.0834 Inf -0.721 0.7848
##
## Results are given on the log (not the response) scale.
## P value adjustment: dunnettx method for 3 tests
## Site rate SE df asymp.LCL asymp.UCL
## TP 31.2 1.86 Inf 27.8 35.1
## PFCA 1 28.5 1.69 Inf 25.4 32.0
## PFCA 2 38.5 1.96 Inf 34.8 42.5
## PFCA 3 29.4 1.71 Inf 26.2 33.0
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## # A tibble: 547 × 4
## # Groups: Site, Transect [46]
## Site Transect SPP6 Number_Seedlings
## <chr> <fct> <chr> <dbl>
## 1 CONTROL 1 DIGISC 6
## 2 CONTROL 1 SETFAB 4
## 3 CONTROL 10 ALOCAR 1
## 4 CONTROL 10 CONCAN 1
## 5 CONTROL 3 ALOCAR 1
## 6 CONTROL 6 CONCAN 1
## 7 CONTROL 6 DIGISC 1
## 8 CONTROL 7 OXADIL 1
## 9 CONTROL 7 PANDIC 1
## 10 CONTROL 8 TAROFF 1
## # ℹ 537 more rows
## # A tibble: 40 × 3
## # Groups: Site [4]
## Site Transect Species_richness
## <fct> <fct> <int>
## 1 TP 1 14
## 2 TP 10 13
## 3 TP 2 11
## 4 TP 3 11
## 5 TP 4 12
## 6 TP 5 10
## 7 TP 6 4
## 8 TP 7 12
## 9 TP 8 9
## 10 TP 9 15
## # ℹ 30 more rows
##
## Call:
## glm(formula = Species_richness ~ Site, family = "poisson", data = Uni_Seedling_Diversity)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.40695 0.09492 25.359 <2e-16 ***
## SitePFCA 1 0.22494 0.12729 1.767 0.0772 .
## SitePFCA 2 0.30775 0.12503 2.461 0.0138 *
## SitePFCA 3 0.20312 0.12791 1.588 0.1123
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 53.060 on 39 degrees of freedom
## Residual deviance: 46.581 on 36 degrees of freedom
## AIC: 230.36
##
## Number of Fisher Scoring iterations: 4
## Analysis of Deviance Table (Type III tests)
##
## Response: Species_richness
## LR Chisq Df Pr(>Chisq)
## Site 6.4783 3 0.09052 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast estimate SE df z.ratio p.value
## PFCA 1 - TP 0.225 0.127 Inf 1.767 0.1919
## PFCA 2 - TP 0.308 0.125 Inf 2.461 0.0382
## PFCA 3 - TP 0.203 0.128 Inf 1.588 0.2674
##
## Results are given on the log (not the response) scale.
## P value adjustment: dunnettx method for 3 tests
## Site rate SE df asymp.LCL asymp.UCL
## TP 11.1 1.05 Inf 9.22 13.4
## PFCA 1 13.9 1.18 Inf 11.77 16.4
## PFCA 2 15.1 1.23 Inf 12.87 17.7
## PFCA 3 13.6 1.17 Inf 11.50 16.1
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## # A tibble: 40 × 3
## # Groups: Site [4]
## Site Transect Species_richness_native
## <fct> <fct> <int>
## 1 TP 1 12
## 2 TP 10 11
## 3 TP 2 10
## 4 TP 3 10
## 5 TP 4 11
## 6 TP 5 9
## 7 TP 6 4
## 8 TP 7 11
## 9 TP 8 8
## 10 TP 9 15
## # ℹ 30 more rows
##
## Call:
## glm(formula = Species_richness_native ~ Site, family = "poisson",
## data = Uni_SB_Diversity_native)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.313e+00 9.950e-02 23.241 <2e-16 ***
## SitePFCA 1 -3.015e-02 1.418e-01 -0.213 0.832
## SitePFCA 2 8.947e-11 1.407e-01 0.000 1.000
## SitePFCA 3 4.832e-02 1.391e-01 0.347 0.728
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 55.512 on 39 degrees of freedom
## Residual deviance: 55.189 on 36 degrees of freedom
## AIC: 227.08
##
## Number of Fisher Scoring iterations: 4
## Analysis of Deviance Table (Type III tests)
##
## Response: Species_richness_native
## LR Chisq Df Pr(>Chisq)
## Site 0.32366 3 0.9555
## contrast estimate SE df z.ratio p.value
## PFCA 1 - TP -0.0302 0.142 Inf -0.213 0.9820
## PFCA 2 - TP 0.0000 0.141 Inf 0.000 1.0000
## PFCA 3 - TP 0.0483 0.139 Inf 0.347 0.9504
##
## Results are given on the log (not the response) scale.
## P value adjustment: dunnettx method for 3 tests
## Site rate SE df asymp.LCL asymp.UCL
## TP 10.1 1.00 Inf 8.31 12.3
## PFCA 1 9.8 0.99 Inf 8.04 11.9
## PFCA 2 10.1 1.00 Inf 8.31 12.3
## PFCA 3 10.6 1.03 Inf 8.76 12.8
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = everything_mat_hell_bray ~ Site * Type, data = everything_mat_labs, method = "bray", p.adjust = "bonferroni")
## Df SumOfSqs R2 F Pr(>F)
## Site 3 10.843 0.31373 26.5607 0.001 ***
## Type 2 5.650 0.16347 20.7590 0.001 ***
## Site:Type 6 3.781 0.10940 4.6309 0.001 ***
## Residual 105 14.288 0.41341
## Total 116 34.561 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## $parent_call
## [1] "everything_mat_hell_bray ~ Site * Type , strata = Null , permutations 999"
##
## $`PFCA 1_vs_PFCA 2`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 1.9875 0.15857 16.8704 0.001 ***
## Type 2 3.4405 0.27449 14.6020 0.001 ***
## Site:Type 2 0.7445 0.05939 3.1596 0.001 ***
## Residual 54 6.3616 0.50755
## Total 59 12.5340 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 1_vs_PFCA 3`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 5.1421 0.30178 36.3215 0.001 ***
## Type 2 2.8375 0.16653 10.0215 0.001 ***
## Site:Type 2 1.4148 0.08303 4.9968 0.001 ***
## Residual 54 7.6448 0.44866
## Total 59 17.0392 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 1_vs_TP`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 6.0682 0.37289 48.6998 0.001 ***
## Type 2 2.3205 0.14259 9.3114 0.001 ***
## Site:Type 2 1.5299 0.09401 6.1391 0.001 ***
## Residual 51 6.3548 0.39050
## Total 56 16.2733 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 2_vs_PFCA 3`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 2.6237 0.16259 17.8593 0.001 ***
## Type 2 4.4234 0.27412 15.0549 0.001 ***
## Site:Type 2 1.1567 0.07168 3.9367 0.001 ***
## Residual 54 7.9331 0.49161
## Total 59 16.1368 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 2_vs_TP`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 3.8134 0.24391 29.2764 0.001 ***
## Type 2 3.6334 0.23240 13.9473 0.001 ***
## Site:Type 2 1.5447 0.09880 5.9296 0.001 ***
## Residual 51 6.6430 0.42489
## Total 56 15.6345 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 3_vs_TP`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 2.1080 0.13796 13.5635 0.001 ***
## Type 2 4.0498 0.26505 13.0290 0.001 ***
## Site:Type 2 1.1957 0.07825 3.8468 0.001 ***
## Residual 51 7.9262 0.51874
## Total 56 15.2797 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"
##
## Call:
## lm(formula = Directionality ~ Site_labels, data = trajectory.stats2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10512 -0.02697 -0.00408 0.04412 0.07062
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.470789 0.008716 54.016 < 2e-16 ***
## Site_labelsYoung -0.034596 0.012014 -2.880 0.00476 **
## Site_labelsMiddle -0.016128 0.012014 -1.342 0.18214
## Site_labelsOld -0.010394 0.012014 -0.865 0.38877
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04529 on 113 degrees of freedom
## Multiple R-squared: 0.07287, Adjusted R-squared: 0.04826
## F-statistic: 2.961 on 3 and 113 DF, p-value: 0.03531
## Anova Table (Type III tests)
##
## Response: Directionality
## Sum Sq Df F value Pr(>F)
## (Intercept) 5.9843 1 2917.7332 < 2e-16 ***
## Site_labels 0.0182 3 2.9606 0.03531 *
## Residuals 0.2318 113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## contrast estimate SE df t.ratio p.value
## Young - Remnant -0.0346 0.012 113 -2.880 0.0135
## Middle - Remnant -0.0161 0.012 113 -1.342 0.4009
## Old - Remnant -0.0104 0.012 113 -0.865 0.6994
##
## P value adjustment: dunnettx method for 3 tests
## # A tibble: 21,177 × 5
## # Groups: Site, Transect [39]
## Site Transect Type SPP6 Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 Aboveground ABUTHE 0
## 2 PFCA 1 1 Seed_Bank ABUTHE 0
## 3 PFCA 1 1 Seed_Rain ABUTHE 0
## 4 PFCA 1 10 Aboveground ABUTHE 0
## 5 PFCA 1 10 Seed_Bank ABUTHE 0
## 6 PFCA 1 10 Seed_Rain ABUTHE 0
## 7 PFCA 1 2 Aboveground ABUTHE 0
## 8 PFCA 1 2 Seed_Bank ABUTHE 0
## 9 PFCA 1 2 Seed_Rain ABUTHE 0
## 10 PFCA 1 3 Aboveground ABUTHE 0
## # ℹ 21,167 more rows
## Warning in as_grob.default(plot): Cannot convert object of class list into a
## grob.